# -*- coding: utf-8 -*-
"""
Created on Tue Dec 11 14:09:55 2012

@author: Maxim
"""
import Curves
from airfoil import Airfoil, AirfoilPolar
import numpy as ny
import matplotlib.pyplot as plt
import os
import shlex
import subprocess as sp
import paths

def flap_geometry(airfoil, flapRatio, lipLen, leRad, gap, overlap, defl,vectRatio1,vectRatio2, vectRatio3):
    flapRatio = 1 - flapRatio
    vectRatio = [vectRatio1,vectRatio2,vectRatio3]
    loLipRatio = 1.5
#    rTe = 0.0041
#    offset = 0.005
    rTe = 0.0001
    offset = 0.0
#    lipTe = 3.33e-3
    lipTe = 0.001
    if -overlap>gap:
        overlap = 0.95*gap

    upNode = ny.zeros([4,2])
    loNode = ny.zeros([4,2])

    upCurve = Curves.xyCurve(airfoil.upPts)
    loCurve = Curves.xyCurve(airfoil.loPts)
    thickness = upCurve(flapRatio)-loCurve(flapRatio)
    leRad = leRad * thickness

    upNode[0,0] = flapRatio+lipLen
    upNode[0,1] = upCurve(upNode[0,0])
    
    tan0 = upCurve.tan(upNode[0,0])
    upNode[1,0] = upNode[0,0] - lipLen*vectRatio[0]
    upNode[1,1] = -(upNode[0,0]-upNode[1,0])*tan0 + upNode[0,1]
    
    loNode[0,0] = flapRatio + leRad*loLipRatio
    loNode[0,1] = loCurve(loNode[0,0])
    
    tan1 = loCurve.tan(loNode[0,0])
    loNode[1,0] = loNode[0,0] - leRad*vectRatio[2]
    loNode[1,1] = -(loNode[0,0]-loNode[1,0])*tan1 + loNode[0,1]
    
    upNode[2,0] = flapRatio
    loNode[2,0] = flapRatio
    
    upNode[3,0] = flapRatio
    loNode[3,0] = flapRatio
    
    upNode[3,1] = loCurve(upNode[3,0]) + leRad*vectRatio[1]
    loNode[3,1] = upNode[3,1]
    
    upNode[2,1] = upNode[3,1] + leRad
    loNode[2,1] = upNode[3,1] - leRad

    curve3 = Curves.BezierCurve(upNode,15)
    curve4 = Curves.BezierCurve(loNode,8)
    curve4e = ny.flipud(curve4)
    curve4e[:,2] = curve4[:,2]
    curve4e = curve4e[1:]
    curve4e[:,2] = curve4e[:,2]+1.0
    curve11 = ny.vstack([curve3,curve4e])
    
    curve11offset = Curves.xytCurveOffset(curve11,-offset)
    curve6 = Curves.xyCurveSplit(upCurve.pts,upNode[0,0])
    curve7 = Curves.xyCurveSplit(loCurve.pts,0,loNode[0,0])
    curve8 = Curves.xyCurveSplit(loCurve.pts,loNode[0,0])
    
    curve9, curve10 = Curves.trimTe(airfoil.upPts,curve11offset,lipTe)
    curve7, curve10 = Curves.splitCurveSimple(curve7,curve10,[0,0])

    curve12, circle1, curve13 = Curves.roundCorner(curve7,curve10,rTe)
    
    mainSec = ny.vstack([ny.flipud(curve9),curve12,circle1,ny.flipud(curve13[:,0:2])])
    flapSec1 = ny.vstack([ny.flipud(curve6)[:-1,:],curve11[:-1,0:2],curve8])
    flapSec0 = flapSec1

    flapSec1 = Curves.rotate2D(flapSec1,[0.7,-.1],defl)
    flapAdd = Curves.rotate2D(curve11,[0.7,-.1],defl)
    
    refPt = mainSec[-1]
    moveX1 = ny.min(flapSec1[:,0]) - refPt[0] + overlap
    flapSec1[:,0] = flapSec1[:,0] - moveX1
    flapAdd[:,0]  = flapAdd[:,0]  - moveX1
    #gap
    dy = Curves.adjustGap(flapAdd,refPt,gap)
    flapSec1[:,1] = flapSec1[:,1] + dy
    
    return mainSec,flapSec1#,flapSec0


def flap_geometry2(x=None,args=1):
    
    if x==None:
        x = ny.array([0.15, 0.3, 0.05, 0.5, 0.1, 0.5, 0.02, 0.02])
    
    flapRatio = 0.7
    defl = 35.0
    airfoil = Airfoil()
    airfoil.read_txt('aiaa_cfd_final_optimum.txt')
    offset = 0.005 # main section flap offset
    lipTe  = 0.0033 # trailing edge thickness
    rTe    = 0.01 # trailing edge radius

    exitLipLen   = x[0]
    leRadPos     = x[1]
    entryLipLen  = x[2]
    v1           = x[3]
    v2           = x[4]
    v3           = x[5]
    gap          = x[6]
    overlap      = x[7]
    
    upNode = ny.zeros([4,2])
    loNode = ny.zeros([4,2])

    upCurve = Curves.xyCurve(airfoil.upPts)
    loCurve = Curves.xyCurve(airfoil.loPts)

    thickness = upCurve(flapRatio)-loCurve(flapRatio) #local thickness
    leRad = v2 * thickness # V2 length
    leRadPos = leRadPos * thickness + loCurve(flapRatio)

    upNode[0,0] = flapRatio + exitLipLen
    upNode[0,1] = upCurve(upNode[0,0])
    
    tan0 = upCurve.tan(upNode[0,0])
    upNode[1,0] = upNode[0,0] - exitLipLen*v1
    upNode[1,1] = -(upNode[0,0]-upNode[1,0])*tan0 + upNode[0,1]
    
    loNode[0,0] = flapRatio + entryLipLen
    loNode[0,1] = loCurve(loNode[0,0])
    
    tan1 = loCurve.tan(loNode[0,0])
    loNode[1,0] = loNode[0,0] - entryLipLen*v3
    loNode[1,1] = loNode[0,1] - (loNode[0,0]-loNode[1,0])*tan1
    
    upNode[2,0] = flapRatio
    loNode[2,0] = flapRatio
    upNode[3,0] = flapRatio
    loNode[3,0] = flapRatio
    
    upNode[3,1] = leRadPos
    loNode[3,1] = leRadPos
    upNode[2,1] = upNode[3,1] + leRad
    loNode[2,1] = loNode[3,1] - leRad
    
    curve3 = Curves.BezierCurve(upNode,10)
    curve4 = Curves.BezierCurve(loNode,10)
    curve4e = ny.flipud(curve4)
    curve4e[:,2] = curve4[:,2]
    curve4e = curve4e[1:]
    curve4e[:,2] = curve4e[:,2]+1.0
    curve11 = ny.vstack([curve3,curve4e])
    
    curve11offset = Curves.xytCurveOffset(curve11,-offset)
    curve6 = Curves.xyCurveSplit(upCurve.pts,upNode[0,0])
    curve7 = Curves.xyCurveSplit(loCurve.pts,0,loNode[0,0])
    curve8 = Curves.xyCurveSplit(loCurve.pts,loNode[0,0])
    
    curve9, curve10 = Curves.trimTe(airfoil.upPts,curve11offset,lipTe)
    curve7, curve10 = Curves.splitCurveSimple(curve7,curve10,[0,0])

    curve12, circle1, curve13 = Curves.roundCorner(curve7,curve10,rTe)
    
    mainSec = ny.vstack([ny.flipud(curve9),curve12,circle1,ny.flipud(curve13[:,0:2])])
    flapSec1 = ny.vstack([ny.flipud(curve6)[:-1,:],curve11[:-1,0:2],curve8])
    #flapSec0 = flapSec1

    flapSec1 = Curves.rotate2D(flapSec1,[0.7,-.1],defl)
    flapAdd = Curves.rotate2D(curve11,[0.7,-.1],defl)
    
    refPt = mainSec[-1]
    moveX1 = ny.min(flapSec1[:,0]) - refPt[0] + overlap
    flapSec1[:,0] = flapSec1[:,0] - moveX1
    flapAdd[:,0]  = flapAdd[:,0]  - moveX1
    #gap
    dy = Curves.adjustGap(flapAdd,refPt,gap)
    flapSec1[:,1] = flapSec1[:,1] + dy

    return mainSec, flapSec1

def flap_objective2(x):
    try:
        Flap  = flap_geometry2(x)
        polar = calc_Jpolar(0.06,3e6,Flap)
        clmax = ny.max(polar.cl)
        return -clmax,[],[]
    except:
        return 100,[],[]


def save_javafoil_igs(x,i):
    flap = flap_geometry2(x)
    path = 'D:\\doe_flap%d.iges'%i
    write_igs(flap,path)

def write_igs(flap,pathtosave):
    path = paths.MyPaths()
    path.setRandPrefix()
    tmpJournal = path.getTmpFile('jfscript')
    tmpAfFile = path.getTmpFile('dat')
    write_txt(flap,tmpAfFile)
    jouFile = open(tmpJournal,'wt')
    jouFile.write('Options.Country(0)\nGeometry.Clear()\n')
    jouFile.write('Geometry.Open(/%s/)\n' %tmpAfFile)
    jouFile.write('Geometry.Save("%s")'%pathtosave)
    jouFile.write('Exit()')
    jouFile.close()
    cmd = ('%s -cp %s -jar %s Script=%s'%(path.java,path.mhclasses,path.javafoil,tmpJournal))
    args=shlex.split(cmd,False,os.name=='posix')
    ps=sp.Popen(args,stdin=sp.PIPE,stderr=sp.PIPE,stdout=sp.PIPE)
    output, errors = ps.communicate()
    ps.stderr.close()
    os.remove(tmpJournal)
    os.remove(tmpAfFile)

def write_txt(Flap, filePath):
    file1 = open(filePath,'wt')
    file1.write('flap\n')
    mainSec, flap = Flap
    for pt in mainSec:
        file1.write('%.6f\t%.6f\n'%(pt[0],pt[1]))
    file1.write('9999.9\t9999.9\n')
    for pt in flap:
        file1.write('%.6f\t%.6f\n'%(pt[0],pt[1]))
    file1.close()

def _read_Jpolar(polarPath):
    polar = AirfoilPolar()
    polarFile = open(polarPath,'rt')
    lines = polarFile.readlines()
    polarFile.close()
    del lines[0:5]
    for line in lines:
        if line.strip()!='':
            segLine = line.split()
            try:
                polar.alpha = ny.append(polar.alpha, float(segLine[0]))
                polar.cl    = ny.append(polar.cl, float(segLine[1]))
                polar.cd    = ny.append(polar.cd, float(segLine[2]))
                polar.cm    = ny.append(polar.cm, float(segLine[3]))
                polar.TU    = ny.append(polar.TU, float(segLine[4]))
                polar.TL    = ny.append(polar.TL, float(segLine[5]))
                polar.SU    = ny.append(polar.SU, float(segLine[6]))
                polar.SL    = ny.append(polar.SL, float(segLine[7]))
                polar.LD    = ny.append(polar.LD, float(segLine[8]))
                polar.AC    = ny.append(polar.AC, float(segLine[9]))
                polar.CP    = ny.append(polar.CP, float(segLine[10]))
            except ValueError:
                print 'Polar processing failed'
    return polar

def calc_Jpolar(Mach, Re, Flap):
    path = paths.MyPaths()
    path.setRandPrefix()
    tmpJournal = path.getTmpFile('jfscript')
    tmpAfFile = path.getTmpFile('dat')
    tmpPolar  = path.getTmpFile('pol')
    write_txt(Flap,tmpAfFile)
    jouFile = open(tmpJournal,'wt')
    jouFile.write('Options.Country(0)\nGeometry.Clear()\n')
    jouFile.write('Geometry.Open(/%s/)\n' %tmpAfFile)
    jouFile.write('Options.MachNumber(%.4f)\n'%Mach)
    jouFile.write('Options.StallModel(0)\n')
    jouFile.write('Options.TransitionModel(1)\n')
    jouFile.write('Options.GroundEffect(0)\n')
    jouFile.write('Options.HeightOverSpan(0.5)\n')
    jouFile.write('Options.AspectRatio(0)\n')
    if os.name=='nt':
        jouFile.write('Polar.Analyze(%.0f;%.0f;%.0f;%.0f;%.0f;%.2f;100;100;0;0)\n'%(Re,Re,0,-10,30,1))    
    else:                
        jouFile.write('Polar.Analyze(%.0f:%.0f:%.0f:%.0f:%.0f:%.2f:100:100:0:0)\n'%(Re,Re,0,-10,30,1))
    jouFile.write('Polar.Save(/%s/)\n'%tmpPolar)
    #jouFile.write('Geometry.Save("D:\\codes\\actools\\pyAC\\actools\\Jfoil\\doe_flap%d.iges")'%i)
    jouFile.write('Exit()')
    jouFile.close()
    cmd = ('%s -cp %s -jar %s Script=%s'%(path.java,path.mhclasses,path.javafoil,tmpJournal))
    args=shlex.split(cmd,False,os.name=='posix')
    ps=sp.Popen(args,stdin=sp.PIPE,stderr=sp.PIPE,stdout=sp.PIPE)
    output, errors = ps.communicate()
    ps.stderr.close() 
    polar = _read_Jpolar(tmpPolar)
    polar.M  = Mach
    polar.Re = Re
    os.remove(tmpJournal)
    os.remove(tmpAfFile)
    os.remove(tmpPolar)
    return polar

def flap_objective(x):
    af = Airfoil()
    af.read_txt('GA37A315mod_reworked.txt')
    flapRatio = 0.3 
#    lipLen    = 0.13 #0.1 - 0.27 
#    leRad     = 0.15 #0.05 - 0.25
#    gap       = 0.02 #0.005;0.02
#    overlap   = -0.02 #-0.02 0.02
    defl      = 35
#    vratio1   = 0.73 #0.5 0.95
#    vratio2   = 1.3 #1.05 1.5
#    vratio3   = 0.55 #0.2 0.9
    lipLen = x[0]
    leRad = x[1]
    gap = x[2]
    overlap = x[3]
    vratio1 = x[4]
    vratio2 = x[5]
    vratio3 = x[6]
    try:
        Flap  = flap_geometry(af,flapRatio,lipLen, leRad, gap, overlap, defl,vratio1,vratio2,vratio3)
        polar = calc_Jpolar(0.06,3e6,Flap)
        clmax = ny.max(polar.cl)
        #print clmax
        return -clmax,[],[]
    except:
        return 100,[],[]

def test_func():
    af = Airfoil()
    af.read_txt('GA37A315mod_reworked.dat')
    print af.up, af.coord
    flapRatio = 0.3 
    lipLen    = 0.1643088
    leRad     = 0.14054948
    gap       = 0.0157
    overlap   = -0.02
    defl      = 35
    vratio1   = 0.77033947
    vratio2   = 1.39164569
    vratio3   = 0.55
    Flap = flap_geometry(af, 0.3, lipLen, leRad,gap,overlap,35,vratio1,vratio2,vratio3)
    mainSec, flap = Flap
    polar = calc_Jpolar(0.06,3e6,Flap)
    
    plt.figure(1)
    plt.hold(True)
    plt.plot(mainSec[:,0],mainSec[:,1])
    plt.plot(flap[:,0],flap[:,1])
    plt.axis([0,1.1,-.5,.5])
    plt.grid(True)
    plt.show()

def plot_results():
#    xBest =[0.18978,0.24883,0.00824,0.02996,0.50002,1.37683,0.20097]
    #xBest =[0.18913399,0.13287908,0.00548619,0.02999233,0.50557791,1.49805884,0.20113115]
    xBest=[0.19417077,0.24936138,0.01010853,0.0299996,0.50077753,1.39159708,0.2021984]
    af = Airfoil()
    af.read_txt('GA37A315mod_reworked.txt')
    FlapOpt = flap_geometry(af, 0.3, xBest[0], xBest[1],xBest[2],xBest[3],35,xBest[4],xBest[5],xBest[6])
    mainSec, flap, flap0 = FlapOpt
    flap1 = Curves.rotate2D(flap,[0.785,-0.135],-35)
    mainSec, flap, flap0 = FlapOpt
   # write_txt(FlapOpt,'KLA100flap_20121213.txt')
    plt.figure(1)
    plt.hold(True)
    plt.plot(mainSec[:,0],mainSec[:,1])
    plt.plot(flap[:,0],flap[:,1])
    plt.plot(flap0[:,0],flap0[:,1])
    #plt.axis([0,1.1,-.5,.5])
    plt.axis('equal')
    plt.grid(True)
    plt.show()
    
class flap2igs():

    def __init__(self,mainSec,flapSec,tol=10.**-5,FREECADPATH='/usr/lib/freecad/lib'):
        import sys
        sys.path.append(FREECADPATH)        
        import FreeCAD
        import Draft
        mainSec=self.unique(mainSec,tol)
        flapSec=self.unique(flapSec,tol)
        Mmin=mainSec[:,0].argmin()
        Fmin=flapSec[:,0].argmin()
        
        mainTop=mainSec[0:Mmin+1,:]
        mainBot=mainSec[Mmin:,:]
        flapTop=flapSec[0:Fmin+1,:]
        flapBot=flapSec[Fmin:,:]

        FreeCAD.newDocument('flap2igs')    
        objs=[]    
#        self.mainObj=self.getSpline(mainSec)
#        self.flapObj=self.getSpline(flapSec)
#        objs.append(self.mainObj)
#        objs.append(self.flapObj)
                
#        self.cutObj=Draft.cut(self.boundaryObj,self.mainObj)
#        self.cutObj=Draft.cut(self.boundaryObj,self.mainObj)
        objs.append(self.getSpline(mainTop))
        objs.append(self.getSpline(mainBot))
        objs.append(self.getSpline(flapTop))
        objs.append(self.getSpline(flapBot))
        objs=self.closeTE(objs,mainSec)
        objs=self.closeTE(objs,flapSec)
        self.objs=objs
#    def setBoundary(self,width=50.):
#        p1=tuple([-width/2,-width/2])
#        p2=tuple([-width/2,+width/2])
#        p3=tuple([+width/2,+width/2])
#        p4=tuple([+width/2,-width/2])     
#        return self.getSpline([p1,p2,p3,p4,p1])
    def writeIGS(self,filePath):
        import Part
        Part.export(self.objs,filePath)
    def unique(self,a,tol=0.0):
        import numpy
        shp=numpy.shape(a)
        rm=[]
        for i in range(shp[0]-1):
            val1=a[i,:]
            val2=a[i+1,:]
            test=(numpy.abs(val1-val2)<tol).all()
            if test==True:
                rm.append(i)
        a=numpy.delete(a,rm,0)
        return a
    def getSpline(self,section):
        import Draft
        import FreeCAD
        points=[]        
        for pts in section:
            points.append(FreeCAD.Vector(pts[0],pts[1],0.0))
        return Draft.makeBSpline(points,closed=False,face=False,support=None)
        #return Draft.makeWire(points,closed=True,face=False,support=None)
    def closeTE(self,objs,section):
        import Draft
        if (section[-1,:]<>section[0,:]).any():
            p1=tuple([section[-1,0],section[-1,1],0.0])
            p2=tuple([section[0,0],section[0,1],0.0])
            objs.append(Draft.makeLine(p1,p2))
        return objs


def run_test1():
    lb = ny.array([0.05, 0.1, 0.05, 0.1, 0.1, 0.1, 0.005, -0.03])
    ub = ny.array([0.20, 0.9, 0.20, 0.9, 0.9, 0.9, 0.030,  0.03])
    x0 = lb+ub/2.0
    print flap_objective2(x0)
    flap, main = flap_geometry2(x0)
    plt.figure(1)
    plt.hold(True)
    plt.axis('equal')
    plt.plot(flap[:,0],flap[:,1],'bo-')
    plt.plot(main[:,0],main[:,1],'ro-')
    plt.show()

def run_objective_for_sens():
    fid = open('flap_variables.txt','rt')
    line = fid.readline().split()
    fid.close()
    x = ny.array([float(xx) for xx in line])
    f,g,h = flap_objective2(x)
    fid = open('flap_objective_value.txt','wt')
    fid.write('%.10f'%f)
    fid.close()

if __name__=="__main__":
    run_objective_for_sens()